tw.cl.lm <- function (data, dv, cntrls, ia, fe, c1, c2, formula = NULL, 
          demeaned = FALSE, ...) 
{
  if (is.null(formula)) {
    if (!is.null(fe) & is.null(ia)) {
      data <- na.omit(data[, c(dv, cntrls, fe, c1, c2)])
      if (demeaned == TRUE) {
        library(dplyr)
        library(magrittr)
        data %<>% group_by(.dots = fe) %>% mutate_at(vars(-one_of(fe, 
                                                                  c1, c2)), funs(as.numeric(as.character(.)))) %>% 
          mutate_if(is.numeric, funs(. - mean(., na.rm = T))) %>% 
          ungroup() %>% as.data.frame()
        formula = paste0(dv, "~", paste(cntrls, collapse = "+"))
      }
      else {
        formula = paste0(dv, "~", paste(cntrls, collapse = "+"), 
                         "+ factor(", fe, ")")
      }
    }
    if (!is.null(fe) & !is.null(ia)) {
      ia.cntrl1 <- stringr::str_extract(ia, "\\*.*$")
      ia.cntrl1 <- stringr::str_remove(ia.cntrl1, "\\*")
      ia.cntrl1 <- trimws(ia.cntrl1)
      ia.cntrl2 <- stringr::str_extract(ia, "^.*\\*")
      ia.cntrl2 <- stringr::str_remove(ia.cntrl2, "\\*")
      ia.cntrl2 <- trimws(ia.cntrl2)
      data <- na.omit(data[, c(dv, cntrls, fe, ia.cntrl1, 
                               ia.cntrl2, c1, c2)])
      if (demeaned == TRUE) {
        library(dplyr)
        library(magrittr)
        data %<>% group_by(.dots = fe) %>% mutate_at(vars(-one_of(ia.cntrl1, 
                                                                  ia.cntrl2, fe, c1, c2)), funs(as.numeric(as.character(.)))) %>% 
          mutate_if(is.numeric, funs(. - mean(., na.rm = T))) %>% 
          ungroup() %>% as.data.frame()
        formula = paste0(dv, "~", paste(cntrls, collapse = "+"), 
                         " + ", ia)
      }
      else {
        formula = paste0(dv, "~", paste(cntrls, collapse = "+"), 
                         "+ factor(", fe, ")", " + ", ia)
      }
    }
    if (is.null(fe) & !is.null(ia)) {
      ia.cntrl1 <- stringr::str_extract(ia, "\\*.*$")
      ia.cntrl1 <- stringr::str_remove(ia.cntrl1, "\\*")
      ia.cntrl1 <- trimws(ia.cntrl1)
      ia.cntrl2 <- stringr::str_extract(ia, "^.*\\*")
      ia.cntrl2 <- stringr::str_remove(ia.cntrl2, "\\*")
      ia.cntrl2 <- trimws(ia.cntrl2)
      data <- na.omit(data[, c(dv, cntrls, fe, ia.cntrl1, 
                               ia.cntrl2, c1, c2)])
      formula = paste0(dv, "~", paste(cntrls, collapse = "+"), 
                       " + ", ia)
    }
    if (is.null(fe) & is.null(ia)) {
      data <- na.omit(data[, c(dv, cntrls, c1, c2)])
      formula = paste0(dv, "~", paste(cntrls, collapse = "+"))
    }
  }
  m <- lm(formula, data)
  if (is.null(c1)) {
    vcv <- multiwayvcov::cluster.vcov(m, cbind(data[, c2]))
  }
  if (is.null(c2)) {
    vcv <- multiwayvcov::cluster.vcov(m, cbind(data[, c1]))
  }
  if (!is.null(c1) & !is.null(c2)) {
    vcv <- multiwayvcov::cluster.vcov(m, cbind(data[, c1], 
                                               data[, c2]))
  }
  se <- lmtest::coeftest(m, vcv)
  newList <- list(lm = m, vcv = vcv, tw.se = se)
  return(newList)
}

ia.plot.df <- function (model, effect, moderator, interaction, type, varcov = "default", 
          minimum = "min", maximum = "max", incr = "default", num_points = 50, 
          conf = 0.95, mean = FALSE, median = FALSE, alph = 80) 
{
  if (type == "binary") {
    minimum = 0
    maximum = 1
    incr = 1
  }
  makeTransparent <- function(someColor, alpha = alph) {
    newColor <- col2rgb(someColor)
    apply(newColor, 2, function(curcoldata) {
      rgb(red = curcoldata[1], green = curcoldata[2], 
          blue = curcoldata[3], alpha = alpha, maxColorValue = 255)
    })
  }
  if (varcov == "default") {
    covMat = vcov(model)
  }
  else {
    covMat = varcov
  }
  mod_frame = model.frame(model)
  beta_1 = model$coefficients[[effect]]
  beta_3 = model$coefficients[[interaction]]
  if (minimum == "min") {
    min_val = min(mod_frame[[moderator]])
  }
  else {
    min_val = minimum
  }
  if (maximum == "max") {
    max_val = max(mod_frame[[moderator]])
  }
  else {
    max_val = maximum
  }
  if (min_val > max_val) {
    stop("Error: Minimum moderator value greater than maximum value.")
  }
  if (incr == "default") {
    increment = (max_val - min_val)/(num_points - 1)
  }
  else {
    increment = incr
  }
  x_2 <- seq(from = min_val, to = max_val, by = increment)
  delta_1 = beta_1 + beta_3 * x_2
  var_1 = covMat[effect, effect] + (x_2^2) * covMat[interaction, 
                                                    interaction] + 2 * x_2 * covMat[effect, interaction]
  se_1 = sqrt(var_1)
  z_score = qnorm(1 - ((1 - conf)/2))
  upper_bound = delta_1 + z_score * se_1
  lower_bound = delta_1 - z_score * se_1
  max_y = max(upper_bound)
  min_y = min(lower_bound)
  data <- as.data.frame(cbind(x_2, delta_1, upper_bound, lower_bound))
  if (type == "binary") {
    data[1, 1] <- moderator
    data[2, 1] <- stringr::str_remove(interaction, paste0(effect, 
                                                          ":"))
  }
  return(data)
}

